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ABSTRACT 

We propose new algorithms for computing triangular de- 
compositions of polynomial systems incrementally. With 
respect to previous works, our improvements are based on 
a weakened notion of a polynomial GCD modulo a regu- 
lar chain, which permits to greatly simplify and optimize 
the sub-algorithms. Extracting common work from similar 
expensive computations is also a key feature of our algo- 
rithms. In our experimental results the implementation of 
our new algorithms, realized with the RegularChains library 
in Maple, outperforms solvers with similar specifications 
by several orders of magnitude on sufficiently difficult prob- 
lems. 

1. INTRODUCTION 

The Characteristic Set Method [21] of Wu has freed Ritt's 
decomposition from polynomial factorization, opening the 
door to a variety of discoveries in polynomial system solving. 
In the past two decades the work of Wu has been extended to 
more powerful decomposition algorithms and applied to dif- 
ferent types of polynomial systems or decompositions: dif- 
ferential systems [2] IIP] , difference systems [9], real para- 
metric systems 22 , primary decomposition [17] . cylindrical 
algebraic decomposition 0j. Today, triangular decomposi- 
tion algorithms provide back-engines for computer algebra 
system front-end solvers, such as Maple's solve command. 

Algorithms computing triangular decompositions of polyno- 
mial systems can be classified in several ways. One can first 
consider the relation between the input system 5* and the 
output triangular systems Si, ■ ■ ■ , S e . From that perspec- 
tive, two types of decomposition are essentially different: 
those for which Si,. . . ,S e encode all the points of the zero 
set S (over the algebraic closure of the coefficient field of S) 
and those for which Si, ■ ■ ■ ,S e represent only the "generic 
zeros" of the irreducible components of S. 

One can also classify triangular decomposition algorithms 
by the algorithmic principles on which they rely. From this 



other angle, two types of algorithms are essentially differ- 
ent: those which proceed by variable elimination, that is, by 
reducing the solving of a system in n unknowns to that of a 
system in n — 1 unknowns and those which proceed incre- 
mentally, that is, by reducing the solving of a system in m 
equations to that of a system in m — 1 equations. 

The Characteristic Set Method and the algorithm in [20] 
belong to the first type in each classification. Kalkbrener's 
algorithm [11] , which is an elimination method solving in the 
sense of the "generic zeros" , has brought efficient techniques, 
based on the concept of a regular chain. Other works [121 
116] on triangular decomposition algorithms focus on incre- 
mental solving. This principle is quite attractive, since it 
allows to control the properties and size of the intermedi- 
ate computed objects. It is used in other areas of polyno- 
mial system solving such as the probabilistic algorithm of 
Lecerf [13] based on lifting fibers and the numerical method 
of Sommese, Verschelde, Wample [18] based on diagonal ho- 
motopy. 

Incremental algorithms for triangular decomposition rely on 
a procedure for computing the intersection of an hypersur- 
face and the quasi-component of a regular chain. Thus, the 
input of this operation can be regarded as well-behaved ge- 
ometrical objects. However, known algorithms, namely the 
one of Lazard [12] and the one of the second author [16] are 
quite involved and difficult to analyze and optimize. 

In this paper, we revisit this intersection operation. Let 
R — k[a;i, . . . , x n \ be the ring of multivariate polynomi- 
als with coefficients in k and ordered variables x = xi < 
■ ■ ■ < x n . Given a polynomial p £ R and a regular chain 
T C k[xi, . . . , x n ], the function call lntersect(p, T, R) returns 
regular chains Ti, . . . , T e C k[si, . . . , x n ] such that we have: 



V{p) n W(T) C W(Ti) U ■ • ■ U W(T e ) C V(p) n W(T). 

(See Section [2] for the notion of a regular chain and related 
concepts and notations.) We exhibit an algorithm for com- 
puting lntersect(p, T, R) which is conceptually simpler and 
practically much more efficient than those of [121 116] . Our 
improvements result mainly from two new ideas. 



Weakened notion of polynomial GCDs modulo regular 
chain. Modern algorithms for triangular decomposition rely 
implicitly or explicitly on a notion of GCD for univariate 
polynomials over an arbitrary commutative ring. A formal 
definition was proposed in [16] (see Definition!]]) and applied 



to residue class rings of the form A = k[x]/sat(T) where 
sat(T) is the saturated ideal of the regular chain T. A mod- 
ular algorithm for computing these GCDs appears in [14| : 
if sat(T) is known to be radical, the performance (both in 
theory and practice) of this algorithm are very satisfactory 
whereas if sat(T) is not radical, the complexity of the algo- 
rithm increases substantially w.r.t. the radical case. In this 
paper, the ring A will be of the form k[x]/^sat(T) while our 
algorithms will not need to compute a basis nor a character- 
istic set of \J sat(T). For the purpose of polynomial system 
solving (when retaining the multiplicities of zeros is not re- 
quired) this weaker notion of a polynomial GCD is clearly 
sufficient. In addition, this leads us to a very simple pro- 
cedure for computing such GCDs, see Theorem [1] To this 
end, we rely on the specialization property of subresultants. 
Appendix IA1 reviews this property and provides corner cases 
for which we could not find a reference in the literature. 

Extracting common work from similar computations. Up 

to technical details, if T consists of a single polynomial t 
whose main variable is the same as p, say v, computing 
lntersect(p, T, R) can be achieved by successively computing 

(si) the resultant r of p and t w.r.t. v, 

(S2) a regular GCD of p and t modulo the squarefree part 
of r. 

Observe that Steps (si) and (S2) reduce essentially to com- 
puting the subresultant chain of p and t w.r.t. v. The algo- 
rithms of Section [3] extend this simple observation for com- 
puting lntersect(p, T, R) with an arbitrary regular chain. In 
broad terms, the intermediate polynomials computed during 
the "elimination phasis" of lntersect(p, T, R) are recycled for 
performing the "extension phasis" at essentially no cost. 

The techniques developed for lntersect(p, T, R) are applied 
to other key sub-algorithms, such as: 

• the regularity test of a polynomial modulo the satu- 
rated of a regular chain, see Section [4] 

• the squarefree part of a regular chain, see Appendix [Bl 

The primary application of the operation Intersect is to ob- 
tain triangular decomposition encoding all the points of the 
zero set of the input system. However, we also derive from 
it in Section [5] an algorithm computing triangular decompo- 
sitions in the sense of Kalkbrener. 

Experimental results. We have implemented the algorithms 
presented in this paper within the Regular Chains library in 
Maple, leading to a new implementation of the Triangularize 
command. In Section [71 we report on various benchmarks. 
This new version of Triangularize outperforms the previous 
ones (based on [15]) by several orders of magnitude on suffi- 
ciently difficult problems. Other Maple commands or pack- 
ages for solving polynomial systems (the WSolve package, 
the Groebner: -Solve command and the Groebner : -Basis 
command for a lexicographical term order) are also outper- 
formed by the implementation of the algorithms presented 



in this paper both in terms of running time and, in the case 
of engines based on Grobner bases, in terms of output size. 

2. REGULAR CHAINS 

We review hereafter the notion of a regular chain and its 
related concepts. Then we state basic properties (Proposi- 
tions [TJ[2l[3l[4l and Corollaries [1] [2| of regular chains, which 
are at the core of the proofs of the algorithms of Section 2] 

Throughout this paper, k is a field, K is the algebraic closure 
of k and k[x] denotes the ring of polynomials over k, with 
ordered variables x = x\ < • • • < i„. Let p £ k[x]. 

Notations for polynomials. If p is not constant, then the 
greatest variable appearing in p is called the main variable 
of p, denoted by mvar(p). Furthermore, the leading coeffi- 
cient, the degree, the leading monomial, the leading term 
and the reductum of p, regarded as a univariate polynomial 
in mvar(p), are called respectively the initial, the main de- 
gree, the rank, the head and the tail of p; they are denoted by 
init(p), mdeg(p), rank(p), head(p) and tail(p) respectively. 
Let q be another polynomial of k[x]. If q is not constant, 
then we denote by prem(p, q) and pquo(p, q) the pseudo- 
remainder and the pseudo-quotient of p by q as univariate 
polynomials in mvar(g). We say that p is less than q and 
write p -< q if either p £ k and q £ k or both are non- 
constant polynomials such that mvar(p) < mvar(q) holds, 
or mvar(p) = mvar(g) and mdeg(p) < mdeg(g) both hold. 
We write p ~ q if neither p -< q nor q -< p hold. 

Notations for polynomial sets. Let F C k[x]. We denote by 
(F) the ideal generated by F in k[x]. For an ideal T C k[x], 
we denote by dim(I) its dimension. A polynomial is regular 
modulo I if it is neither zero, nor a zerodivisor modulo X. 
Denote by V(F) the zero set (or algebraic variety) of F in 
K™. Let h £ k[x]. The saturated ideal of I w.r.t. h, denoted 
by X : ft 00 , is the ideal {q £ k[x] 3m £ N s.t. h m q £ I}. 

Triangular set. Let T C k[x] be a triangular set, that is, a 
set of non-constant polynomials with pairwise distinct main 
variables. The set of main variables and the set of ranks of 
the polynomials in T are denoted by mvar(T) and rank(T), 
respectively. A variable in x is called algebraic w.r.t. T 
if it belongs to mvar(T), otherwise it is said free w.r.t. T. 
For v £ mvar(T), denote by T v the polynomial in T with 
main variable v. For v £ x, we denote by T <v (resp. T>„) 
the set of polynomials t £ T such that mvar(t) < v (resp. 
mvar(f) > v) holds. Let hr be the product of the initials 
of the polynomials in T. We denote by sat(T) the saturated 
ideal of T defined as follows: if T is empty then sat(T) is 
the trivial ideal (0), otherwise it is the ideal (T) : ■ The 
quasi-component W(T) of T is defined as V(T) \ V{hr)- 
Denote W(T) = V(sat(T)) as the Zariski closure of W(T). 
For F C k[x], we write Z(F,T) := V{F) n W(T). 

Rank of a triangular set. Let S C k[x] be another tri- 
angular set. We say that T has smaller rank than S and 
we write T -< S if there exists v £ mvar(T) such that 
rank(T<„) = rank(5 , < t) ) holds and: (i) either v £ mvar(S'); 
(ii) or v £ mvar(S') and T v -< S v . We write T ~ S if 
rank(T) = rank(S). 



Iterated resultant. Let p, q G k[x]. Assume g is noncon- 
stant and let v = mvar(g). We define res(p, q, v) as follows: if 
the degree deg(p, v) of p in v is null, then res(p, q, v) = p; oth- 
erwise ies(p, q, v) is the resultant of p and q w.r.t. v. Let T 
be a triangular set of k[x]. We define res(p, T) by induction: 
if T — 0, then res(p, T) — p; otherwise let v be greatest vari- 
able appearing in T, then res(p, T) = res(res(p, T v , v), T< v ). 

Regular chain. A triangular set T C k[x] is a regular chain 
if: (i) either T is empty; (ii) or T\{T max } is a regular chain, 
where T max is the polynomial in T with maximum rank, and 
the initial of T max is regular w.r.t. sat(T \ {T max }). The 
empty regular chain is simply denoted by 0. 

Triangular decomposition. Let F C k[x] be finite. Let 
X := {Ti , . . . , T e } be a finite set of regular chains of k[x] . We 
call X a Kalkbrener triangular decomposition of V(F) if we 
have V(F) = U| =1 W(T 4 ). We call X a Lazanf- triangular 
decomposition of V(F) if we have V(F) = Uf =1 W / (T l ). 

Proposition 1 (Th. 6.1. in |T]). Lei p and T be re- 
spectively a polynomial and a regular chain o/k[x]. Then, 
prem(p,T) = holds if and only if p G sat(T) holds. 

Proposition 2 (Prop. 5 in [16]). LetT andT' betwo 
regular chains of k[x] such that J sat(T) C J sat(T') and 
dim (sat(T)) = dim (sat(T')) hold. Let p g k[x] suc/i that p 
is regular w.r.t. sat(T). Then p is also regular w.r.t. sat(T'). 

Proposition 3 (Prop. 4.4 in [L). iei p g k[x] and 
T C k[x] be a regular chain. Let v — mvar(jp) and r — 
prem(p,T> v ) such that r g w sat(T <v ) holds. Then, we have 
p G y/sat(T). 

Corollary 1. Let T and T' be two regular chains of 
\s[xx, . . . , Xk\, where 1 < k < n. Letp G k[x] with mvar(p) = 
Xk+i such that init(p) is regular w.r.t. both sat(T) and 
sat(T'). Assume that sat(T) C y/sat(T') holds. Then 
we also have \J sat(T U p) C \J sat(T' U p). 

Proposition 4 (Lemma 4 in [3]). Let p e k[x]. Let 
T C k[x] be a regular chain. Then the following statements 
are equivalent: 

(i) the polynomial p is regular w.r.t. sat(T), 
(ii) for each prime ideal p associated with satiT), we have 

(Hi) the iterated resultant res(p, T) is not zero. 

Corollary 2. Let p e k[x] and T C k[x] be a regular 
chain. Let v := mvar(p) and r := res(p,T> v ). We have: 

(1) the polynomial p is regular w.r.t. sat(T) if and only if 
r is regular w.r.t. sa£(T<„); 

(2) if v ^ mvar(T) and init(p) is regular w.r.t. sat(T), 
then p is regular w.r.t. sat(T). 



3. REGULAR GCDS 

As mentioned before, Definition [1] was introduced in [16] 
as part of a formal framework for algorithms manipulating 
regular chains [7] 1121 [5] 1111 123] . In the present paper, the 
ring A will always be of the form k[x]/-\/sat(T). Thus, a 
regular CCD of p, t in A[y] is also called a regular CCD of 
p, t modulo \J sat(T). 

Definition 1. Let A be a commutative ring with unity. 
Let p,t,g G A[y] with t =^ and g ^ 0. We say that g £ A[y] 
is a regular GCD of p, t if: 

(Ri) the leading coefficient of g in y is a regular element; 

(R2) g belongs to the ideal generated by p and t in A[y]; 

(R3) if deg(g,y) > 0, then g pseudo-divides both p and t, 
that is, prem(p, g) = prem(t, g) — 0. 

Proposition 5. For 1 < k < n, let T c k[asi, . . . , £ fc _i] 
be a regular chain, possibly empty. Let p,t,g g k[a?i, . . . , Xk] 
be polynomials with main variable x^. Assume T U {t} is 
a regular chain and g is a regular GCD of p and t modulo 
y/sat(T). We have: 

(i) if mdeg(g) = mdeg(t), then \J sat(T U t) = \J sat(T U g) 
andW(TUt) C Z(h g , T U t) U W(T U g) both hold, 

(ii) if mdeg(g) < mdeg(t), let q = pquo(t, g), then T U q is 
a regular chain and the following two relations hold: 

(ii.a) y/sat(T U t) = sat(T U g) n sj ' sat(T U q), 
(li.b) W(TUt) C Z(h g ,TUt) U W(TUg)L)W(TUq), 

(Hi) W(TUg) C V(p), 

(iv) Z(p,TUt) C W(TU g) U Z({p,h g },TUt). 

Proof. We first establish a relation between p, t and g. 
By definition of pseudo-division, there exist polynomials q, r 
and a nonnegtive integer eo such that 

h g °t = qg + r and r g Vsat(T) (1) 

both hold. Hence, there exists an integer ei > such that: 

(h T r(h e g H~q 9 r G (T) (2) 

holds, which implies: t G ^/sat(T U g). We first prove (i). 
Since mdeg(i) = mdeg((/) holds, we have g G k[si, . . . , x k -i], 
and thus we have h g ° ht = qh g . Since h t and h g are reg- 
ular modulo sat(T'), the same property holds for q. To- 
gether with ([2]), we obtain g G -y/sat(T U t). Therefore 
\J sat(T U t) = \J sat(T U g). The inclusion relation in (i) 
follows from ([1]). 

We prove (ii). Assume mdeg(t) > mdeg(g). With (JTJ 
and (J2| , this hypothesis implies that T U q is a regular chain 
and t G ^/sat(T U g) holds. Since t G ^sat(T U g) also 
holds, v /sat(T U t) is contained in \/sat(TU.g) n v /sat(T U g). 
Conversely, for any / G \l sat(T U g) PI \J sat(T U g), there 
exists an integer e-i > and a G k[x] such that (h g h q ) e2 f 2 — 



aqg G sat(T) holds. With Q we deduce that / G v /sat(TUt) 
holds and so does (ii.a). With we have (ii.b) holds. 

We prove (Hi) and (iv). Definition [T] implies: prem(p, g) G 
N /sat(T). Thus p G \/sat(T U g) holds, that is, W(TUg) C 
V(p), which implies (in). Moreover, since g G (p, i, y sat(T)), 
we have Z(p, T U t) C V^ff), so we deduce (iu). □ 

Let p, i be two polynomials of k[asi, . . . , Xk], for k > 1. Let 
m = deg(p,a;fc), n = mdeg(t, x^). Assume that m,n > 1. 
Let A = min(m, n). Let T be a regular chain of k[a?i, . . . , Kfc-i]. 
Let B = k[xi, . . . , Xk-i] and A = B/ x /sat(T). 

Let So, • ■ • , S\-\ be the subresulant polynomials [151 [8] of 
p and t w.r.t. Xk in B[sfe]. Let = coeff(5'i, x\) be the 
principle subresultant coefficient of Si, for < i < A — 1. 
If m > n, we define 5a = t, Sa+i = p, sa = init(f) and 
sa+i = init(p). If m < n, we define S\ = p, Sa+i = t, 
s\ — init(p) and sa+i = init(t). 

The following theorem provides sufficient conditions for Sj 
(with 1 < j < A + 1) to be a regular GCD of p and t in 
A[x k ]. 

Theorem 1. Let j be an integer, with 1 < j < A + 1, 
such that Sj is a regular element of A and such that for any 
< i < j, we have Si = in A. Then Sj is a regular GCD 
of p and t in A[xk]. 

Proof. By Definition [1] it suffices to prove that both 
prem(p, Sj , Xk) = and prem(f, Sj, Xk) = hold in A. By 
symmetry we only prove the former equality. 

Let p be any prime ideal associated with sat(T). Define 
D = k[a;i, . . . ,Xk-i]/p and let L be the fraction field of the 
integral domain D. Let (f> be the homomorphism from B 
to L. By Theorem |3] of Appendix \K\ we know that 4>(Sj) 
is a GCD of (f>{p) and 4>(t) in L[xfe]. Therefore there ex- 
ists a polynomial q of h[xk] such that p = qSj in L[a;fe], 
which implies that there exists a nonzero element a of D 
and a polynomial q' of B[xfe] such that ap = q' Sj in B[a;fe]. 
Therefore prem(ap, Sj) — in D[xft], which implies that 
prem(p, Sj) = in D[j:fc]. Therefore prem(p, Sj) belongs to 
p and thus to N /sat(T). So prem(p, Sj, x k ) = in A. □ 

4. THE INCREMENTAL ALGORITHM 

In this section, we present an algorithm to compute Lazard- 
Wu triangular decompositions in an incremental manner. 
We recall the concepts of a process and a regular (delayed) 
split, which were introduced as Definitions 9 and 11 in [16| . 
To serve our purpose, we modify the definitions as below. 

Definition 2. A process of k[x] is a pair (p,T), where 
p G k[x] is a polynomial and T C k[x] is a regular chain. 
The process (0, T) is also written as T for short. Given 
two processes (p,T) and (p',T r ), let v and v' be respectively 
the greatest variable appearing in (p,T) and (p',T'). We 
say (p,T) -( (p',T') if: (i) either v < v' ; (ii) or v — v' 
and dimT < dimT'; (iii) or v = v' , dimT = dimT' and 
T -< T'; (iv) or v = v' , dimT = dimT', T ~ T' and 



p<p . We write (p,T) ~ (p',T') if neither (p,T) X (p',T') 
nor (p',T r ) -< (p,T) hold. Clearly any sequence of processes 
which is strictly decreasing w.r.t. -< is finite. 

Definition 3. Let Ti, 1 < i < e, be regular chains of 
k[x]. Let p G k[x]. We call Ti,...,T e a regular split of 
(p, T) whenever we have 



(Li) ^Jsat{T) C ^Jsat{T t ) 

(L 2 ) W(Ti) C V(p) (or equivalents p G ^Jsat{Ti)) 
(L 3 ) V(p)nW(T)CUt =1 W(Ti) 

We write as (p, T) — )• Ti, . . . ,T e . Observe that the above 
three conditions are equivalent to the following relation. 

V(p) n W(T) c W{Ti) u • ■ • u W(r e ) c v( P ) n W(T). 

Geometrically, this means that we may compute a little more 
thanV(p)nW{T); however, W(Ti)U- ■ -UW(T e ) is a "sharp" 
approximation of the intersection ofV(p) and W(T). 

Next we list the specifications of our triangular decompo- 
sition algorithm and its subroutines. We denote by R the 
polynomial ring k[x], where x = a;i < • • • < x n . 

Triangularize(i ? , R) 

• Input: F, a finite set of polynomials of R 

• Output: A Lazard-Wu triangular decomposition of 
V(F). 

Intersect(p, T, R) 

• Input: p, a polynomial of R; T, a regular chain of R 

• Output: a set of regular chains {Ti, . . . , T e } such chat 
(p,T) — >Ti,...,T B . 

Regularize^, T, R) 

• Input: p, a polynomial of 7?; T, a regular chain of R. 

• Output: a set of pairs {[pi, Ti], . . . , [p e , T e ]} such that 
for each i,l < i < e: (1) Ti is a regular chain; (2) 
p = Pi mod yJsaX(Ti)\ (3) if pi = 0, then p t G v /sat(T 4 ) 
otherwise p t is regular modulo N /sat(T i ); moreover we 
have T — > T\, . . . ,T e . 

SubresultantChain(p, q, v, R) 

• Input: v, a variable of {xi, . . . , x„}; p and q, polyno- 
mials of R, whose main variables are both v. 

• Output: a list of polynomials (So, ■ ■ ■ , S\), where 
A = min(mdeg(p), mdeg(g)), such that Si is the i-th 
subresultant of p and q w.r.t. v. 
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Input size 


Output size 






#v 


#e 


dcg 


dim 


GL 


GS 


GD 


TL 


TK 


1 


4corps-lparameter-homog 


4 


3 


8 


1 


- 


- 


21863 


- 


30738 


2 


8-3-config-Li 


12 


7 


2 


7 


67965 


- 


72698 


7538 


1384 


3 


Alonso-Li 


7 


4 


4 


3 


1270 


- 


614 


2050 


374 


4 


Bezicr 


5 


3 


6 


2 


_ 


_ 


32054 


_ 


114109 


5 


Cheaters- homotopy- 1 


7 


3 


7 


4 


26387452 


- 


17297 


- 


285 


7 


childDraw-2 


10 


10 


2 





938846 


_ 


157765 


_ 


_ 


8 


Cinquin-Demongeot-3-3 


4 


3 


4 


1 


1652062 


_ 


680 


2065 


895 


9 


Cinquin-Demongeot-3-4 


4 


3 


5 


1 


_ 


_ 


690 


_ 


2322 


10 


collins-jsc02 


5 


4 


3 


1 


_ 


_ 


28720 


2770 


1290 


11 


f-744 


12 


12 


3 


1 


102082 


- 


83559 


4509 


4510 


12 


Haas5 


4 


2 


10 


2 


- 


- 


28 


- 


548 


14 


Lichtblau 


3 


2 


11 


1 


6600095 


_ 


224647 


110332 


5243 


16 


Liu-Lorcnz 


5 


4 


2 


1 


47688 


123965 


712 


2339 


938 


17 


Mehta2 


11 


8 


3 


3 


- 


- 


1374931 


5347 


5097 


18 


Mehta3 


13 


10 


3 


3 


- 


- 


- 


25951 


25537 


19 


Mehta4 


15 


12 


3 


3 


- 


- 


- 


71675 


71239 


21 


p3p- isosceles 


7 


3 


3 


4 


56701 




1453 


9253 


840 


22 


p3p 


8 


3 


3 


5 


160567 




1768 




1712 


23 


Pavelle 


8 


4 


2 


4 


17990 




1552 


3351 


1086 


24 


Solotareff-4b 


5 


4 


3 


1 


2903124 




14810 


2438 


872 


25 


Wang93 


5 


4 


3 


1 


2772 


56383 


1377 


1016 


391 


26 


Xia 


6 


3 


4 


3 


63083 


2711 


672 


1647 


441 


27 


xy-5-7-2 


6 


3 


3 


3 


12750 




599 




3267 



Table 1 The input and output sizes of systems 



RegularGcd(p, q, v, S, X, R) 



CleanChain(C, X, Xi, R) 



• Input: v, a variable of {xi, . . . ,x n }, 

— X, a regular chain of R such that mvar(T) < v, 

— p and q, polynomials of R with the same main 
variable v such that: init(g) is regular modulo 
^/sat(X); res(p, q, v) belongs to -^/sat(X), 

— S, the subresultant chain of p and q w.r.t. v. 

• Output: a set of pairs {[<?i, Xi], . . . , [g e , X e ]} such that 
X — > Xi, . . . , X e and for each X;: if dimX = dimX» , 
then gi is a regular GCD of p and q modulo sat(Xi); 
otherwise gi = 0, which means undefined. 

lntersectFree(p, Xi, C, R) 



• Input: Xi, a variable of x; p, a polynomial of R with 
main variable Xi\ C, a regular chain of k[a;i, . . . , Xj— x]- 

• Output: a set of regular chains {Xi, . . . , X e } such that 
(p,C) — ► (Xi,...,X e ). 



lntersectAlgebraic(p, X, Xi, S, C, R) 

• Input: p, a polynomial of R with main variable Xi, 

— X, a regular chain of R, where Xi € mvar(X), 

— S, the subresultant chain of p and T Xi w.r.t. Xi, 

— C, a regular chain of k[xi, . . . , Xi-i], such that: 
inh^Xc J is regular modulo ^/sat(C); the resul- 
tant of p and T Xi , which is So, belongs to \J sat(C). 

• Output: a set of regular chains Xl, . . . ,T e such that 
(p,CUT Xi ) — >Xi,...,X e . 



Input: X, a regular chain of R; C, a regular chain of 
k[a;i, . . . , Xi-i] such that y/ sa.t(T <Xi ) C y' sat(C). 

Output: if Xi £ mvar(X), return C; otherwise return 
a set of regular chains {Xi, . . . , X e } such that inh^X^) 
is regular modulo each sat(Xj), ^/sat(C) C -^/sat(Xj) 
and W{C) \ V{init(T Xi )) C U| =1 W(Tj). 



Extend(C,X,a;i,i?) 



Input: C, is a regular chain of k[xi, 



X, a 



regular chain of R such that y / sat(T <Xi ) C ^/sat(C). 

Output: a set of regular chains {Xi,...,X e } of R 
such that W(CUT> Xi ) C U| =1 W(Tj) and ^/sat(T) C 



Algorithm SubresultantChain is standard, see [8]. The algo- 
rithm Triangularize is a principle algorithm which was first 
presented in [16j . We use the following conventions in our 
pseudo-code: the keyword return yields a result and ter- 
minates the current function call while the keyword output 
yields a result and keeps executing the current function call. 

5. PROOF OF THE ALGORITHMS 

Theorem 2. All the algorithms in Fig. 1 terminate. 

Proof. The key observation is that the flow graph of Fig. 
1 can be transformed into an equivalent flow graph satisfying 
the following properties: (1) the algorithms Intersect and 
Regularize only call each other or themselves; (2) all the other 
algorithms only call either Intersect or Regularize. Therefore, 
it suffices to show that Intersect and Regularize terminate. 



Algorithm 1: lntersect(p, T, R) 

1 if prem(p, T) = then return {T}; 

2 if p G k then return { }; 

3 r:=p;P:= {r}; S := { }; 

4 while mvarir) € mvariT) do 

5 ii := mvar(r); src := SubresultantChain(r, T v , v, R); 

6 S := S U {src}; r := resultant(src); 

7 if r = then break; 

8 if r £ k then return { }; 

9 _ P := P U {r} 

10 X:={0}; X':={};i:=l; 
while i < n do 
for C 6 X do 

if Xi mvar(P) and Xi ^ mvar(T) then 
| X' := X'UCIeanChain(C,T,a: i+ i,.R) 
else if £i ^ mvar(P) then 
| X' := X' U CleanChain(CUT Xi ,T,a; i+ i,fl) 
else if Xi rnvar(T) then 

for D G \ntersectFree(P Xi , Xi, C, R) do 
|_ X' := X' U CleanChain(L>,T,;r l+ i,F) 

else 

for D e lntersectAlgebraic(P It ,T,x iy S Xi , C, R) do 
|_ %' := X' U CleanChain(D,T, x l+ i,R) 



Algorithm 5: Regularize^, T, R) 



ll 

12 
13 
14 
15 
16 
17 
18 
19 

20 
21 
22 



23 



X:=X'; X' := { }; i :=i + l 



24 return X 



Algorithm 3: lntersectFree(p, Xi, C, R) 



Algorithm 4: lntersectAlgebraic(p, T, Xi, S, C, R) 

for [g,D] G Regu\arGcd(p,T Xt ,Xi,S,C,R) do 
if dim D < dim C then 

for E G CleanChain(D, T,Xi,R) do 
|_ output lntersectAlgebraic(p, T, Xi, S, E, R) 

else 

output D U g; 

for E G lntersect(mii(gi), D, R) do 
for F G OeanChain (E,T,x i, R) do 
|_ output lntersectAlgebraic(p, T, x t , S, F, R) 



if p G k or T = then return [p, T]; 

w := mvar(p); 

if v (f: mvar(T) then 

for [f,C] G Regularize(imi(p), T,R) do 

if / = then output Regularize(tail(p), C, 7?); 



else output [p, C]; 



8 else 



Algorithm 2: RegularGcd(p, q, v, S, T, R) 

1 X:={(T,1)}; 

2 while X / do 

3 let (C*,i) G X; X := X\ {(C,i)}; 

4 for [/, D] G Regularize^, C, R) do 

5 if dimD < dimC then output [0, D] ; 

6 else if / = then X := X U {(£>, z + 1)} ; 

7 else output [Si , D] 



9 
10 
11 
12 
13 

14 
15 
16 
17 
18 
19 

20 
21 

22 
23 
24 
25 
26 
27 



src := SubresultantChain(p, T v ,v,R); r := resultant(src); 
for [f,C] G Regularize^, T <V ,R) do 
if dimC < dimT<v then 

for D G Extend (C, T, v, R) do 
|_ output Regularize^, D, R) 

else if / / then output [p, C U T> v ] ; 
else 

for [g,D] G RegularGcd(p, T v , v, src, C, R) do 
if dim D < dimC then 

for E G Extend{D,T,v,R) do 
|_ output Regularize^, E, R); 

else 

if mdeg(p) = mdeg(Ft,) then output 
[0,DUT>„]; next; 
output [0,flUjU T>„]; 
Q := pquo(T„,.g); 

output Regularize^, D U q U T>„, F); 
for F G lntersect(/i 9 , F/, R) do 
for F G Extend(F,T,u,F) do 
j_ output Regularize^, F, R) 



Algorithm 6: Extend (C, T, xt, R) 



for [/, D] G Regularize(imi(p), C, R) do 

if / = then output lntersect(tail(p), D, R) ; 
else 

output D U p; 

for E G lntersect(mii(p), _D, R) do 
[_ output Intersect (tail (p), E, R) 



1 if T> :Ei = then return C; 

2 let p G T with greatest main variable; T' := T \ {p}; 

3 for D G Extend(C,T',2; l ,i?) do 



for [/, E] G Regularize(imt(p), D) do 
|_ if / / then output E U p; 



Algorithm 7: CleanChain(C, T, Zi, 7?) 



1 if ^ mvar(F) or dim C = dim T <Xi then return C; 

2 for [f,D] G Regularized^ J, C,R) do 

3 |_ if / 7^ then output D 



Algorithm 8: TriangularizefF, i?) 
l if F = { } then return {£ 



2 Choose a polynomial p G F with maximal rank; 

3 for T G Triangularize(F \ {p}, R) do 

4 |_ output lntersect(p, T, R) 



Intersect 




Regularize IntersectAlgebraic 

/ u\ / 

Extend RegularGcd 

Li 

Figure 1: Flow graph of the Algorithms 

Note that the input of both functions is a process, say (p, T). 
One can check that, while executing a call with (p,T) as 
input, any subsequent call to either functions Intersect or 
Regularize will take a process (p',T') as input such that 
(p', T') -< (p, T) holds. Since a descending chain of processes 
is necessarily finite, both algorithms terminate. □ 

Since all algorithms terminate, and following the flow graph 
of Fig. 1, each call to one of our algorithms unfold to a finite 
dynamic acyclic graph (DAG) where each vertex is a call to 
one of our algorithms. Therefore, proving the correctness of 
these algorithms reduces to prove the following two points. 

• Base: each algorithm call, which makes no subsequent 
calls to another algorithm or to itself, is correct. 

• Induction: each algorithm call, which makes subse- 
quent calls to another algorithm or to itself, is correct, 
as soon as all subsequent calls are themselves correct. 

For all algorithms in Fig. 1, proving the base cases is straight- 
forward. Hence we focus on the induction steps. 

Proposition 6. IntersectFree satisfies its specification. 
Proof. We have the following two key observations: 

• C — > Di, . . . ,D S , where Di are the regular chains in 
the output of Regularize. 

• V(p) n W{D) = W{D,p) U V(init(p),tail(p)) n W(D). 
Then it is not hard to conclude that (p, C) — > Ti, . . . , T e . □ 

Proposition 7. IntersectAlgebraic is correct. 

PROOF. We need to prove: (p,C UT Xi ) — 5- Ti, . . . , T e . 
Let us prove (Li) now, that is, for each regular chain Tj in 
the output, we have -\/sat(C U T Xi ) C *J sat(Tj). First by 
the specifications of the called functions, we have -y/sat(C) C 
^/sat(£>) C x /sat(_B), thus, sJs&t(C U T Xi ) C v /sat(_B U T Xq ) 
by Corollary [1] since init(T Xi ) is regular modulo both sat(C) 



and sat(_E). Secondly, since g is a regular GCD of p and T Xi 
modulo y/ sat(D), we have -\/sat(C U T Xi ) C y' sat(D U g) 
by Corollaries [T] and Proposition [5] 

Next we prove (L2 ). It is enough to prove that W(D U 
g) != V(p) holds. Since g is a regular GCD of p and T Xi 
modulo ^Jsa£[D), the conclusion follows from point (Hi) of 
Proposition [S] 

Finally we prove (L 3 ), that is Z{p, C U T Xi ) C \J° W(Tj). 
Let Di , . . . , D s be the regular chains returned from Algo- 
rithm RegularGcd. We have C — s- D\, . . . , D a , which im- 
plies Z(p, C U T Xi ) C Uj = iZ(p, Dj U T Xi ). Next since g is a 
regular GCD of p and T Xi modulo \Js&t(Dj), the conclusion 
follows from point (iv) of Proposition [5] □ 

Proposition 8. Intersect satisfies its specification. 

Proof. The first while loop can be seen as a projection 
process. We claim that it produces a nonempty triangular 
set P such that V(p) n W(T) = V(P) n W{T). The claim 
holds before staring the while loop. For each iteration, let P' 
be the set of polynomials obtained at the previous iteration. 
We then compute a polynomial r, which is the resultant of 
a polynomial in P' and a polynomial in T. So r £ (P',T). 
By induction, we have (p,T) = (P,T). So the claim holds. 

Next, we claim that the elements in X satisfy the follow- 
ing invariants: at the beginning of the i-th iteration of the 
second while loop, we have 

(1) each C £ X is a regular chain; if T Xi exists, then 
init(T Xi ) is regular modulo sat(C), 

(2) for each C G X, we have \/s&t(T <Xi ) C x /sat(C), 

(3) for each C £ X, we have W{C) C V(P <Xi ), 

(4) V(p)nW(T)C[J 0e% Z(P> Xi ,CUT> Xi ). 

When i = n+1, we then have v /sat(T) C v /sat(C), W(C) C 
V{P) C V(p) for each C £ X and V(p)rW(T) C UoeiW(C). 
So (Lr), (£2), (£3) of Definition [3] all hold. This concludes 
the correctness of the algorithm. 

Now we prove the above claims (1), (2), (3), (4) by induc- 
tion. The claims clearly hold when i = 1 since C — 
and V(p) n W(T) = V(P) n W{T). Now assume that the 
loop invariants hold at the beginning of the i-th iteration. 
We need to prove that it still holds at the beginning of the 
(i + l)-th iteration. Let C £ X be an element picked up at 
the beginning of i-th iteration and let L be the set of the 
new elements of X' generated from C. 

Then for any C £ L, claim (1) clearly holds by specification 
of CleanChain. Next we prove (2). 

• if Xi ^ mvar(T), then T <Xi+1 = T <Xi . By induction 
and specifications of called functions, we have 

^Js&t{T <Xi+1 ) C v / sat(C) C ^s&t(C). 



• if Xi S mvar(T), by induction we have *J sa,t(T <Xi ) C 
yj sat(C) and init(T Xi ) is regular modulo both sat(C) 
and sat(T<a, i ). By Corollary [T] we have 

^sat{T <Xi+1 ) C V , sat(Cur li ) C ^sat{C). 

Therefore (2) holds. Next we prove claim (3). By induction 
and the specifications of called functions, we have W(C) C 
W(CUT X i ) C Vj P <a!< ). Secondly, we have W^J C ^(PsJ. 
Therefore W(C) C V(P <a!i+1 ), that is (3) holds. Finally, 
since V(P^)nW/(CUT^)\V(init(T^ +1 )) C U c , ejL VF(C"), 
we have Z(P> Xi ,Cl>T> Xi ) C U , gi Z(P> ai+1 , C' UT> ai+1 ), 
which implies that (4) holds. This completes the proof. □ 

Proposition 9. Regularize satisfies its specification. 

Proof. If v ^ mvar(T), the conclusion follows directly 
from point (2) of Corollary [2] From now on, assume v G 
mvar(T). Let L be the set of pairs [p',7 1 '] in the output. We 
aim to prove the following facts 

(1) each T' is a regular chain, 

(2) if p = 0, then p is zero modulo \J sat(T'), otherwise p 
is regular modulo sat(T), 

(3) we have v /sat(T) C x /sat(T'), 

(4) we have W(T) C U T , gL VF(T'). 

Statement (1) is due to Proposition [2] Next we prove (2). 
First, when there are recursive calls, the conclusion is obvi- 
ous. Let [/, C] be a pair in the output of Regularize^, T <V ,R). 
If / 7^ 0, the conclusion follows directly from point (1) of 
Corollary [2] Otherwise, let [g, D] be a pair in the output of 
the algorithm Regu larGcd (p, T v , v, src, C, R). If mdeg(g) = 
mdeg(T„), then by the algorithm of RegularGcd, g = T v . 
Therefore we have prem(p, T v ) G ^/sat(C), which implies 
that p G \Js&t(C U T> v ) by Proposition [3] 

Next we prove (3). Whenever Extend is called, (3) holds 
immediately. Otherwise, let [/, C] be a pair returned by 
Regularize^, T <v , R). When / 7^ 0, since v/sat(T< v ) C 
\J sat(C) holds, we conclude w sat(T) C ^/sat(C U T>„) 
by Corollary [1] Let [g,D] G RegularGcd(p, T v , v, src, C, R). 
Corollary[T]and point (ii) of Proposition0imply that sat(T) 
N /sat(D ur>„), x/sat(r) C x /sat(D U g U T>„) together 
with x /sat(T) C ^sat(D U g U T>„) hold. Hence (3) holds. 

Finally by point (ii.6) of Proposition[5] we have W(DUT V ) C 
Z(ft 9 ,DUT„)UW(DUj)Ulf(DU9). So (4) holds. □ 

Proposition 10. Extend satisfies its specification. 

Proof. It clearly holds when T> Xi = 0, which is the 
base case. By in duction and the specification of Regularize, 
we know that \J sat(T') C J sat(_B). Since init(p) is reg- 
ular modulo both sat(T') and sat(P), by Corollary [T] we 



have \J sat(T) C sj sat(_B U p). On the other hand, we have 
W(CUT^ X .) C UW(_D) and \ V(h p ) C U 

Therefore #(C U T> a4 ) C U] =1 W(T f ), where Tr, . . . , T e are 
the regular chains in the output. □ 

Proposition 11. CleanChain satisfies its specification. 
Proof. It follows directly from Proposition [2] □ 

Proposition 12. RegularGcd satisfies its specification. 

Proof. Let [gi,Ti\, i = l,...,e, be the output. First 
from the specification of Regularize, we have T — > T\, . . . ,T e . 
When dimTi = dimT, by Proposition and Theorem [T] gi 
is a regular CCD of p and q modulo \J sat(T). □ 

6. KALKBRENER DECOMPOSITION 

In this section, we adapt the Algorithm Triangularize (Algo- 
rithm |8}, in order to compute efficiently a Kalkbrener trian- 
gular decomposition. The basic technique we rely on follows 
from Krull's principle ideal theorem. 

Theorem 3. LetF C k[x] be finite, with cardinality #(P) . 
Assume F generates a proper ideal o/k[x]. Then, for any 
minimal prime ideal p associated with (F) , the height of p is 
less than or equal to #(F). 

Corollary 3. Let 1 be a Kalkbrener triangular decom- 
position ofV(F). Let T be a regular chain ofl, the height 
of which is greater than #(F). Then T\ {T} is also a Kalk- 
brener triangular decomposition ofV(F). 

Based on this corollary, we prune the decomposition tree 
generated during the computation of a Lazard-Wu triangu- 
lar decomposition and remove the computation branches in 
which the height of every generated regular chain is greater 
than the number of polynomials in F. 

Next we explain how to implement this tree pruning tech- 
nique to the algorithms of Section[4] Inside Triangularize, de- 
fine A = #(P) and pass it to every call to Intersect in order 
to signal Intersect to output only regular chains with height 
no greater than A. Next, in the second while loop of Inter- 
sect, for the i-th iteration, we pass the height A — #(T> Xi+1 ) 
to CleanChain, IntersectFree and IntersectAlgebraic. 

In IntersectFree, we pass its input height A to every function 
call. Besides, Lines 5 to 6 are executed only if the height of 
D is strictly less than A, since otherwise we would obtain 
regular chains of height greater than A. In other algorithms, 
we apply similar strategies as in Intersect and IntersectFree. 

7. EXPERIMENTATION 

Part of the algorithms presented in this paper are imple- 
mented in Maple14 while all of them are present in the cur- 
rent development version of Maple. Tables 1 and 2 report 
on our comparison between Triangularize and other Maple 
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Table 2 Timings of Triangularize versus other solvers 



solvers. The notations used in these tables are defined be- 
low. 

Notation for Triangularize. We denote by TK and TL the 
latest implementation of Triangularize for computing, respec- 
tively, Kalkbrener and Lazard-Wu decompositions, in the 
current version of Maple. Denote by TK14 and TL14 the 
corresponding implementation in Maple14. Denote by TK13, 
TL13 the implementation based on the algorithm of |16j in 
MAPLE13. Finally, STK and STL are versions of TK and 
TL respectively, enforcing that all computed regular chains 
are squarefree, by means of the algorithms in Appendix [B] 

Notation for the other solvers. Denote by GL, GS, GD, re- 
spectively the function Groebner:-Basis (plex order), Groebner:- 
Solve, Groebner:-Basis (tdeg order) in current beta version of 
Maple. Denote by WS the function wsolve of the package 
Wsolve [19j ■ which decomposes a variety as a union of quasi- 
components of Wu Characteristic Sets. 

The tests were launched on a machine with Intel Core 2 
Quad CPU (2.40GHz) and 3.0Gb total memory. The time- 
out is set as 3600 seconds. The memory usage is limited to 
60% of total memory. In both Table 1 and 2, the symbol 
"-" means either time or memory exceeds the limit we set. 

The examples are mainly in positive dimension since other 
triangular decomposition algorithms are specialized to di- 
mension zero [6]. All examples are in characteristic zero. 

In Table 1, we provide characteristics of the input systems 
and the sizes of the output obtained by different solvers. 
For each polynomial system F C Q[x], the number of vari- 
ables appearing in F, the number of polynomials in F, the 
maximum total degree of a polynomial in F, the dimension 
of the algebraic variety V(F) are denoted respectively by 
#e, deg, dim. For each solver, the size of its output is 
measured by the total number of characters in the output. 



To be precise, let "dec" and "gb" be respectively the out- 
put of the Triangularize and Groebner functions. The Maple 
command we use are length(convert(map(Equations, dec, R), 
string)) and length(convert(gb, string)). From Table 1, it is 
clear that Triangularize produces much smaller output than 
commands based on Grobner basis computations. 

TK, TL, GS, WS (and, to some extent, GL) can all be seen 
as polynomial system solvers in the sense of that they pro- 
vide equidimensional decompositions where components are 
represented by triangular sets. Moreover, they are imple- 
mented in Maple (with the support of efficient C code in 
the case of GS and GL) . The specification of TK are close to 
those of GS while TL is related to WS, though the triangular 
sets returned by WS are not necessarily regular chains. 

In Table 2, we provide the timings of different versions of 
Triangularize and other solvers. From this table, it is clear 
that the implementations of Triangularize, based on the al- 
gorithms presented in this paper (that is TK14, TL14, TK, 
TL) outperform the previous versions (TK13, TL13), based 
on [16] . by several orders of magnitude. We observe also that 
TK outperforms GS and GL while TL outperforms WS. 
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APPENDIX 

A. SPECIALIZATION PROPERTIES OF SUB- 
RESULTANT CHAINS 

Let A be a commutative ring with identity and let k < I 
be two positive integers. Let M be an k x I matrix with 
coefficients in A. Let Mi be the square submatrix of M 
consisting of the first k— 1 columns of M and the i t h column 
of M, for i = k ■ ■ ■ I. Let det(Mi) be the determinant of 
Mi. We denote by dpol(M) the element of A[x], called the 
determinant polynomial of M, given by 



det M k x + det M k +ix 



detM e . 



Let fi(x), . . . , fk{x) be a set of polynomials of A[x]. Let 
I = 1 + max(deg fi(x), . . . , deg fk(x)). The matrix M of 
/i, . . . , fk is defined by Mij = coeff(/ M x e ~ 3 ). 

Let / = a m x m + ■ ■ ■ + do, g = b„x" + • • • + bo be two 
polynomials of A[x] with positive degrees m and n. Let 
A = min(m, n). For any < i < A, let M be the matrix of 
the polynomials x n ~ 1 ~ l f, . . . , xf, f, x m ~ 1 ~ % g, . . . , xg, g. We 
define the i t h subresultant of / and g, denoted by Si(f,g), 
as 



Si(f,g) 



dpolOr™- 1 
dpol(M). 



X J, /, X 



■,xg,g) 



Note that Si(f,g) is a polynomial in A[x] with degree at 
most i. Let Si(f,g) — coeS(Si(f, g), x r ) and call it the prin- 
ciple subresultant coefficient of Si. 

Let B be a UFD. Let 4> be a homomorphism from A to B, 
which induces naturally also a homomorphism from A[x] to 
M[x]. Let m' = deg(</>(/)) and n' = deg(0(g)). 



Lemma 1. For any integer < k < A, if (f)(sk) / 0, 
then <f>{a m ) and <f>(b n ) does not vanish at the same time. 
Moreover, we have both deg^(/) > k and deg 0(g) > k. 



Proof. Observe that 



Sk = 



do 



Om ■ ■ ■ ftfc 

b n b n -l ■ ■ ■ bo 



bn b n -l 



bk 



Therefore there exists i > k, j > k such that 4( a i) 7^ and 
4>(bj) 7^ 0. The conclusion follows. □ 



Done. □ 



Lemma 2. Assume that 4( s o) 
Then, if m < n, we have 



(sa-i) = 0. 



(1) if 4>{a m ) 7^ and 4(b n ) = ■ ■ ■ = 4{bm) = 0, then 4(g) = 


(2) if4(a m ) = and 4(b n ) ± 0, then </>(f) = 
Symmetrically, if m > n, we have 

(3) if(f>(b n ) 7^ and 4( a m) = ■ ■ • = <f>(a n ) — 0, then <j){f) = 


(4) if 0(fe„) = and 0(a m ) / 0, then 4(g) = 

Proof. We prove (1) and (2), whose correctness implies 
(3) and (4) by symmetry. Let i = A ~ 1 = m — 1, then we 
have 



5 m _i = dpol(a;"-"7, . . . , xf, f,g). 



Therefore 



Sm — 1 — 



(Ml 



bm bm—1 



So from 4>{b n ) = • • • = 4>{b m ) = and 0(s m _i) = 0, we 
conclude that <j>(b m -i) = 0. On the other hand, if 4( a m) = 
and 4>(b n ) 7^ 0, then 4( a m-i) = 0. 



Now let consider S m -2- We have 



Sm-2 = 



(Ml 



b n 



bn 



Am a m — l 
bm-1 bm-2 

bm-1 bm-2 



From 4>{b m -\) — 0, we conclude that 4(bm-2) = 0. From 
4>(a m -i) = 0, we conclude that 4( a m-2) = 0. 

So on so forth, finally, if 4( a m) 7^ and 0(6„) = ■ ■ ■ = 
4>{bm) = 0, we deduce that 4(bi) = 0, for all < i < m — 1, 
which implies that 0(g) = 0; if 0(a m ) = and 4(bn) 7^ 0, 
we deduce that 4( a m-i) — • • • — 4( a o) = 0, which implies 
that 4(f) = 0. □ 

Lemma 3. Let i be an integer such that 1 < i < A. As- 
sume that (j}(a m ) 7^ 0. If i < n , t/ien uie /ia?;e 

= 0(a m )™-™'dpol(x"'- 1 -V(/), . . . ,x<f>(f), <f>(f), 

Proof. If i < n', then n — n' < n — i. Therefore we have 

4(Si) = 4>(&pcl(x n - 1 - i f,...,xf,f,x m - 1 - i g,...,xg,g)) 
= ^(dpoKa;"- 1 - - • - , z<H/), 

= 0(a m )"-"'dpol( 2 ;"'- 1 -V(/), ■ • ■ , atf (/), <H/), 
z™- 1 -^),...,^),^))) 



Theorem 4. W^e /ia«e f/ie following relations between the 
subresultants and the GCD of 4>(f) and cj>(g): 

1. Let k, < k < X, be an integer such that 4( s k) 7^ and 
for any i, < i < k, 4( s i) — 0. Then gcd(0(/), 4(g)) = 

4(s k ). 

2. Assume that cj>( s i) = for all < i < A. we have the 
following cases 

(a) if m < n and 4( a m) 7^ 0, then gcd(4(f),4(g)) = 
4(f); symmetrically, if m > n and 4(bn) 7^ 0, then 
we have gcd(0(/), 4(g)) = 4(g) 

(b) if m < n and 4( a m) = but 0(6„) 7^ 0, then we 
have gcd(4(f) , 4(g)) ~ 4(g); symmetrically, ifm > 
n and 4(bn) ~ but 4( a m) 7^ 0, then we have 

gcd(4(f),4(g)) = 4(f) 

(c) if 4(a m ) = 4(bn) = 0, then 

gcd(4(f),4(g)) = gcd(4(red(f)),4(red(g))) 

Proof. Let us first prove (1). W.l.o.g, we assume 4( a ™) 7^ 
0. From Lemma[T] we know that k < n' . Therefore for i < k, 
we have i < n' . By Lemma [3] 

= 4(am) n - n 'dpo\(x- n '- 1 - i 4(f),...,x4(f),4(f), 

x m - 1 - i 4(g),...,x4(g),4(g)) 

If i < n', we have 0(S S ) = 4M n ~ n ' S,(4(f), 4(g))- H 
i = n', since i < m, we have 

0(5<) = 0(a m )™-™'dpol(x m ~ 1 -V(ff),-.-,^(ff)>(3)) 
= 0(a m )-">(V) m - 1 -V((?). 

So for all i < k, we have Si(4(f) , 4(g)) = 0. If < n', we 
have s fc (0(/), ^(ff)) / 0. So gcd(<M/), <^(ff)) = 4(S k ). If ft = 
n', we have 4{°n') = 4(bk) 0. Therefore gcd(0(/), 0(<7)) = 

0( 3 ) = 4(S k ). 

Next we prove (2a). By symmetry, we prove it when m < n. 
If 4(bn) = • ■ ■ = 4(bm) = 0, it follows directy from Lemma[2] 
Otherwise, we have n' > m. By Lemma O for alH < m we 
have 

0(5 4 ) = 4(am) n - n 'dpol(x n '- 1 - i 4(f),...,x4(f),4(f), 
x m ~ 1 -'4(g),..., x4(g), 4(g)) 

That is 4(Si) = 4M n ~ n 'S l (4(f), 4(g))- Since 4( Si ) = 0, 
we deduce that 4(Si) = gcd(0(/), 4(d))- 

Finally (26) follows directly from Lemma [2] and (2c) is ob- 
viouly true. All done. □ 



B. SQUAREFREE DECOMPOSITION 

Throughout this section, we assume that the coefficient field 
k is of characteristic zero. We propose two strategies for 
computing a squarefree triangular decomposition. The first 
one is a post-processing which applies Algorithm [TT] to every 
regular chain returned by Algorithm [8] The second consists 
of ensuring that, each output or intermediate regular chain 
generared during the execution of Algorithm [8] is squarefree. 



To implement the second strategy, we add an squarefree op- 
tion to Algorithm [8] and each of its subalgorithms. If the 
option is set to true, this option requires that each output 
regular chain is squarefree. This is achieved by using Algo- 
rithm [5] whenever we need to construct new regular chains 
from a previous regular chain T and a polynomial p such 
that T U p is known to be a regular chain. 

Algorithm 9: Squarefree(p, Xi, T, R) 

Input: a polynomial ring R = k[xi, . . . , x n ], a variable 
Xi of R, a squarefree regular chain T of k[a;i, . . . , a 
polynomial p of R with main variable Xi such that T U p is 
a regular chain. 

Output: a set of squarefree regular chains Ti, . . . , T e such 

thatpUT — >T 1: ...,T e . 
p := SquarefreePart(p); 
if mdeg(p) = 1 then return T U p; 
else 

src := SubresultantChain(p,p', Xi, R); 
return Squarefree(p, Xi, src, T, R); 
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Algorithm 10: Squarefree(p, x%, src, T, R) 

Input: a polynomial ring R — k[a;i, . . . , x n ], a variable 
Xi of R, a squarefree regular chain T of k[xi, . . . , Xi-i], a 
squarefree polynomial p of R with main variable Xi such 
that T U p is a regular chain, the sub-resultant chain src of 
p and p' w.r.t Xj. 

Output: a set of squarefree regular chains Ti, . . . , T e such 

thatpUT — >Ti,...,T e . 
r := resultant(src, R); 

i :=U; 

for C G Regularize^, T, i?) do 

if r ^ sat(C) then output C U p; next; 
else 

if dimC = dimT then 
[_Z:= TU {C}; next; 

else 

for [/,£>] G Regularize(imi(p), C,i?) do 
|_ if / / then % :=ZU {D}; 



while T { } do 

let C G T; T := T\ {C}; 

for [g,-D] G Regu larGcd (p, p', Xj, src, C, R) do 
if dim D = dim C then 
output D U pquo(p, g); 
for E G lntersect(inii(g), D, 7?) do 

for [f,F] G Regularize(inii(p), £7, _R) do 
|_ if / / then T := X U {F}; 

else 

for [/, E] G Regularize(imt(p), _D, i?) do 
L if / / then T := 2 U {£}; 



Algorithm 11: Squarefree(T, R) 



Input: a polynomial ring R = k[xi, . . . , x„], a regular 
chain T of R. 



that T — >Ti,...,T e . 

1 T :— {SquarefreePart(p) | p G T}; 

2 5:={}; 

3 for p G T do 



if mdeg(p) > 1 then 
I S :— S U {SubresultantChain(p,p', mvar(p), i?)}; 



6 X:= {0}; X' := { }; i := 1; 

7 while i < n do 



8 
9 
10 
11 
12 
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for C G % do 

if Xi ^ mvar(T) then 
| X' :=T'uCleanChain(C*,T,a; 4+ i,i?) 
else 

if mdeg(T Xi ) = 1 then 

| X' := X'uCleanChain(C , U{T^J,T,a; l+ i, J R) 
else 

for D G Squarefree(T Xi , Xi,S Xi ,C,R) do 
[_ X' := X' U CleanChain(D,T,a; l+ i,_R) 



X:=X'; X' := { }; i := i + 1; 



18 return X 



